function Mat= vdm22(d_col,d_row)
% d_row represent the collotation points degree;
% d_col represent the degree of bform basis
m_row = (d_row + 1)*(d_row + 2)/2;
m_col = (d_col + 1)*(d_col + 2)/2;
[I_row,J_row,K_row] = indices(d_row);
[I_col,J_col,K_col] = indices(d_col);
IM_row = diag(I_row)*ones(m_row,m_col);
JM_row = diag(J_row)*ones(m_row,m_col);
KM_row = diag(K_row)*ones(m_row,m_col);
IM_col = diag(I_col)*ones(m_col,m_row);
JM_col = diag(J_col)*ones(m_col,m_row);
KM_col = diag(K_col)*ones(m_col,m_row);
Mat = (IM_row/d_row).^(IM_col').*(JM_row/d_row).^(JM_col').*(KM_row/d_row).^(KM_col');
IF = gamma(I_col+1);
JF = gamma(J_col+1);
KF = gamma(K_col+1);
A = factorial(d_col)*ones(m_row,m_col)*diag(1./(IF.*JF.*KF));
Mat = A.*Mat;